# using Monte Carlo simulation to simulate the moment of inertia.

from vpython import *
N = 1000
n = 0
R = 1
M = 1
dm = M/N
I = 0

while n < N:
    rtemp = vector(R * (2 * random() - 1), R * (2 * random() - 1), 0)
    if mag(rtemp) < R:
        sphere(pos = rtemp, radius = R/30)
        I = I + dm * mag(rtemp) ** 2
        n = n + 1

print("I = ", I)
print("I_theory = ", 0.5 * M * R ** 2)